recipes %>% 
  ggplot(aes(x=as.factor(rating), fill=as.factor(rating) )) +  
  geom_bar( ) +
  scale_fill_manual(values = c("red4", "red3", "orangered", "orange", "gold", "greenyellow", "green3", "green4") ) +
  theme(legend.position="none") +
  scale_y_continuous(breaks=seq(0,10000,1000)) +
  labs(x = "Rating", y = "Count", 
       title = "Overview of recipes' ratings")

# As we can see, most of the ratings have value above 3.75, more in particular we notice that most of the recipes are rates with 4.375.
# A zero rating indicates that the recipe has not been evaluated yet.

recipes_boxplot <-
  recipes %>% 
  pivot_longer(cols = c(calories, protein, fat, sodium),
               names_to = "nutrition",
               values_to = "n_value")

recipes_calories <- recipes_boxplot %>% 
  filter(nutrition == "calories")

recipes_calories %>%
  ggplot(aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  ggtitle("Boxplot of calories nutritional value") +
  xlab("") +
  ylab("Value")

# We notice that there are recipes with more than 30'000'000 calories which are clearly outliers. We must then discard those values in order to continue with a meaningful analysis. Otherwise from a visual point of view we could not extract any relevant information. If we take out these recipes, we still have roughly 30 recipes which have more than 7000 calories. For the purpose of this precise visualisation, we want to exclude them to have a better overview of the boxplot, but we don't consider those as outliers.

recipes_calories %>% 
  select(title, nutrition, n_value) %>% 
  arrange(desc(n_value)) 

recipes_calories <- recipes_calories %>% 
  filter(n_value <= 7000)


# Calories boxplot

recipes_calories %>%
  ggplot(aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  scale_y_continuous(breaks=seq(0,7000,500)) +
  ggtitle("Boxplot of calories nutritional value") +
  xlab("") +
  ylab("Calories") 

# Calories histogram

recipes_calories %>% 
  ggplot(aes(x=n_value)) +
  geom_histogram(binwidth=50, fill="red3", color="red3", alpha=0.9) +
  theme(plot.title = element_text(size=15)) +
  scale_x_continuous(breaks=seq(0,10000,1000)) +
  scale_y_continuous(breaks=seq(0,1750,250)) +
  ggtitle("Distribution of calories across all recipes") +
  xlab("Calories") +
  ylab("Count") 

recipes_protein <- recipes_boxplot %>% 
  filter(nutrition == "protein")

recipes_protein %>%
  ggplot(aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  ggtitle("Boxplot of protein nutritional value") +
  xlab("") +
  ylab("Value")

# We notice that there are recipes with more than 50'000 grams of protein which are clearly outliers. We must then discard those values in order to continue with a meaningful analysis. Otherwise from a visual point of view we could not extract any relevant information. By checking on the epicurious website recipes  with protein values above 1000, we also verified that the amount of proteins was not justified. We came to that conclusion by evaluating the average values of protein per 100grams of each ingredient in the specific recipe.

recipes_protein %>% 
  select(title, nutrition, n_value) %>% 
  arrange(desc(n_value)) 

recipes_protein <- recipes_protein %>% 
  filter(n_value <= 1000)


# Proteins boxplot

recipes_protein %>%
  ggplot( aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  scale_y_continuous(breaks=seq(0,7000,25)) +
  ggtitle("Boxplot of protein nutritional value") +
  xlab("") +
  ylab("Proteins") 

# Proteins histogram

recipes_protein %>% 
  ggplot(aes(x=n_value)) +
  geom_histogram(binwidth=7, fill="red3", color="red3", alpha=0.9) +
  theme(plot.title = element_text(size=15)) +
  scale_x_continuous(breaks=seq(0,1000,25)) +
  scale_y_continuous(breaks=seq(0,7000,250)) +
  ggtitle("Distribution of proteins across all recipes") +
  xlab("Proteins") +
  ylab("Count") 

recipes_sodium <- recipes_boxplot %>% 
  filter(nutrition == "sodium")

recipes_sodium %>%
  ggplot(aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  ggtitle("Boxplot of sodium nutritional value") +
  xlab("") +
  ylab("Value")

# We notice that there are recipes with more than 100'000 milligrams of sodium which are clearly outliers. We must then discard those values in order to continue with a meaningful analysis. By conducting further research, we realize that sodium values above 30'000 are highly suspicious.

recipes_sodium %>% 
  select(title, nutrition, n_value) %>% 
  arrange(desc(n_value)) 

recipes_sodium <- recipes_sodium %>% 
  filter(n_value <= 30000)


# Sodium boxplot

recipes_sodium %>%
  ggplot(aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  scale_y_continuous(breaks=seq(0,30000,500)) +
  ggtitle("Boxplot of sodium nutritional value") +
  xlab("") +
  ylab("Sodium") 

# Sodium histogram

recipes_sodium %>% 
  ggplot(aes(x=n_value)) +
  geom_histogram(binwidth=50, fill="red3", color="red3", alpha=0.9) +
  theme(plot.title = element_text(size=15)) +
  scale_x_continuous(breaks=seq(0,30000,1000)) +
  scale_y_continuous(breaks=seq(0,1750,250)) +
  ggtitle("Distribution of sodium across all recipes") +
  xlab("Sodium") +
  ylab("Count") 

recipes_fat <- recipes_boxplot %>% 
  filter(nutrition == "fat")

recipes_fat %>%
  ggplot(aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  ggtitle("Boxplot of fat nutritional value") +
  xlab("") +
  ylab("Value")

# We notice that there are recipes with more than 44'000 grams of fat which are clearly outliers. We must then discard those values in order to continue with a meaningful analysis. By checking on the epicurious website recipes  with fat values above 1000, we also verified that the amount of proteins was not justified. We came to that conclusion by evaluating the average values of protein per 100grams of each ingredient in the specific recipe.

recipes_fat %>% 
  select(title, nutrition, n_value) %>% 
  arrange(desc(n_value)) 

recipes_fat <- recipes_fat %>% 
  filter(n_value <= 40000)


# Fat boxplot

recipes_fat %>%
  ggplot( aes(x=nutrition, y=n_value, fill=nutrition)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_light() +
    theme(legend.position="none",
          plot.title = element_text(size=11)) +
  scale_y_continuous(breaks=seq(0,3000,100)) +
  ggtitle("Boxplot of fat nutritional value") +
  xlab("") +
  ylab("Fat") 

# Fat histogram

recipes_fat %>% 
  ggplot(aes(x=n_value)) +
  geom_histogram(binwidth=7, fill="red3", color="red3", alpha=0.9) +
  theme(plot.title = element_text(size=15)) +
  ggtitle("Distribution of fat across all recipes") +
  scale_x_continuous(breaks=seq(0,3000,100)) +
  scale_y_continuous(breaks=seq(0,7000,250)) +
  xlab("Fat") +
  ylab("Count") 

recipes <- recipes %>% 
  filter(calories <= 7000, protein <= 1000, sodium <= 30000, fat <= 40000)


# Scatterplot of Rating-Calories
recipes %>% 
  ggplot(aes(x=calories, y=rating)) +
  geom_point() +
  ggtitle("Scatterplot of rating against calories") +
  xlab("Calories") +
  ylab("Rating") 

# We can observe that the recipes with more than 2000 calories tend to have a higher rating. For instance, few recipes with less than a 3 star rating have more than 2000 calories.


# Scatterplot of Rating-Protein
recipes %>% 
  ggplot(aes(x=protein, y=rating)) +
  geom_point() +
  ggtitle("Scatterplot of rating against proteins") +
  xlab("Proteins") +
  ylab("Rating") 

# We can observe that the recipes with more than 125 grams of proteins tend to have a higher rating. For instance, few recipes with less than a 3 star rating have more than 125 grams of proteins.


# Scatterplot of Rating-Fat
recipes %>% 
  ggplot(aes(x=fat, y=rating)) +
  geom_point() +
  ggtitle("Scatterplot of rating against fat") +
  xlab("Fat") +
  ylab("Rating") 

# We can observe that the recipes with more than 100 grams of fat tend to have a higher rating. For instance, few recipes with less than a 3 star rating have more than 100 grams of fat.


# Scatterplot of Rating-Sodium
recipes %>% 
  ggplot(aes(x=sodium, y=rating)) +
  geom_point() +
  ggtitle("Scatterplot of rating against sodium") +
  xlab("Sodium") +
  ylab("Rating") 

# We can observe that the recipes with more than 5000 milligrams of sodium tend to have a higher rating. For instance, few recipes with less than a 3 star rating have more than 5000 mg of sodium.

corr_nutritional_values = recipes %>% 
  select(rating, calories, protein, fat, sodium) %>% 
  cor()

corrplot(corr_nutritional_values)

# We decide to plot together the variables which highlight a great level of correlation.

# Grouped scatter of calories and fat
recipes_plot1 <- recipes %>% 
  filter(fat <= 400, calories <= 6000)

recipes_plot1 %>% 
  ggplot(aes(x=calories, y=fat, color=rating)) +
  geom_point() + 
  scale_color_gradientn(colours = rainbow(5))

# Grouped scatter of calories and protein
recipes_plot2 <- recipes %>% 
  filter(protein <= 500, calories <= 6000)

recipes_plot2 %>% 
  ggplot(aes(x=calories, y=protein, color=rating)) +
  geom_point() + 
  scale_color_gradientn(colours = rainbow(5))

# Grouped scatter of protein and fat
recipes_plot3 <- recipes %>% 
  filter(fat <= 400, protein <= 350)

recipes_plot3 %>% 
  ggplot(aes(x=protein, y=fat, color=rating)) +
  geom_point() + 
  scale_color_gradientn(colours = rainbow(5))

# Grouped scatter of protein and sodium
recipes_plot4 <- recipes %>% 
  filter(sodium <= 400, protein <= 350)

recipes_plot4 %>% 
  ggplot(aes(x=protein, y=sodium, color=rating)) +
  geom_point() + 
  scale_color_gradientn(colours = rainbow(5))

ingredients_df_bin <- ingredients_df %>%
  mutate(vegetables_bin = as.numeric(if_any(all_of(vegetables_vec), ~.x == 1, na.rm = TRUE)),
         meats_bin = as.numeric(if_any(all_of(all_meats), ~.x == 1, na.rm = TRUE)),
         fish_bin = as.numeric(if_any(all_of(fish_vec), ~.x == 1, na.rm = TRUE)),
         seafood_bin = as.numeric(if_any(all_of(seafood_vec), ~.x == 1, na.rm = TRUE)),
         herbs_bin = as.numeric(if_any(all_of(herbs_vec), ~.x == 1, na.rm = TRUE)),
         nuts_bin = as.numeric(if_any(all_of(nuts_vec), ~.x == 1, na.rm = TRUE)),
         fruits_bin = as.numeric(if_any(all_of(fruits_vec), ~.x == 1, na.rm = TRUE)),
         #drinks_bin = as.numeric(if_any(all_of(drink_no_alcohol_vec), ~.x == 1, na.rm = TRUE)),
         cheese_bin = as.numeric(if_any(all_of(cheeses_vec), ~.x == 1, na.rm = TRUE)),
         dairy_bin = as.numeric(if_any(all_of(dairy_vec), ~.x == 1, na.rm = TRUE)),
         spices_bin = as.numeric(if_any(all_of(spices_vec), ~.x == 1, na.rm = TRUE)),
         cereals_bin = as.numeric(if_any(all_of(cereals_vec), ~.x == 1, na.rm = TRUE)),
         carbs_bin = as.numeric(if_any(all_of(carbs_vec), ~.x == 1, na.rm = TRUE)),
         dessert_bin = as.numeric(if_any(all_of(dessert_vec), ~.x == 1, na.rm = TRUE)),
         egg_bin = (egg)
         ) %>% 
  select(ID, title, contains("bin"), everything())

####testing if I still need to include the general category to create the binary column now that I modified the df to only include recipes with ingredients specified
# 
# #6586
# ingredients_df %>%
#   mutate(vegetables_bin = as.factor(as.numeric(if_any(c(vegetable, all_of(vegetables_vec)), ~.x == 1, na.rm = TRUE)))) %>%
#   filter(vegetables_bin == 1)
# 
# #6586
# ingredients_df %>%
#   mutate(vegetables_bin = as.factor(as.numeric(if_any(all_of(vegetables_vec), ~.x == 1, na.rm = TRUE)))) %>%
#   filter(vegetables_bin == 1)

ingredients_df_total <- ingredients_df %>%
  mutate(total_ingredients = rowSums(select(., c(all_of(all_ingredients)))),
         total_vegetables = rowSums(select(., c(all_of(vegetables_vec)))),
         total_meat = rowSums(select(., c(all_of(all_meats)))),
         total_fish = rowSums(select(., c(all_of(fish_vec)))),
         total_seafood = rowSums(select(., c(all_of(seafood_vec)))),
         total_herbs = rowSums(select(., c(all_of(herbs_vec)))),
         total_nuts = rowSums(select(., c(all_of(nuts_vec)))),
         total_fruits = rowSums(select(., c(all_of(fruits_vec)))),
         #total_drinks = rowSums(select(., c(all_of(drink_no_alcohol_vec)))),
         total_cheese = rowSums(select(., c(all_of(cheeses_vec)))),
         total_dairy= rowSums(select(., c(all_of(dairy_vec)))),
         total_spices= rowSums(select(., c(all_of(spices_vec)))),
         total_cereals= rowSums(select(., c(all_of(cereals_vec)))),
         total_carbs = rowSums(select(., c(all_of(carbs_vec)))),
         total_dessert = rowSums(select(., c(all_of(dessert_vec))))
         ) %>% 
  select(ID, title, contains("total"), everything())

total_join <- ingredients_df_total %>% 
  select(ID, contains("total"))
  
ingredients_df_full <- ingredients_df_bin %>% 
  left_join(total_join) %>% 
  select(ID, title, rating, contains("bin"), contains("total"), everything())

#creating a vector with colnames of all the binary columns to be able to select them more easily afterwards
binary_columns <- ingredients_df_bin %>% 
  select(contains("bin")) %>% 
  colnames()

#adding binary columns to ingredients_df
total_categories <- ingredients_df_bin %>% 
  select(ID, all_of(binary_columns)) %>% 
  pivot_longer(-ID, names_to = "category", values_to = "binary_value") %>% 
  group_by(category) %>% 
  summarise(total = sum(binary_value))

#plotting the frequency of binary columns
total_categories %>%
  ggplot(aes(x=reorder(category,total), y=total, fill=total)) +
  geom_bar(stat = "identity") +
  scale_x_discrete(guide = guide_axis(n.dodge=3))+
  scale_fill_viridis() +
  labs(x = "Category", y = "Amount of recipes", title = "Total amount of recipes containing at least one ingredient in defined categories")

ingredients_df_bin %>% 
  select(contains("bin")) %>%
  mutate(across(everything(), as.factor)) %>% 
  plot_bar(order_bar = FALSE)

ingredients_df_bin %>% 
  select(contains("bin"), rating) %>% 
  plot_correlation()

#testing if corr matrix changes when we set the binary columns as factor instead of num
ingredients_df_bin %>% 
  select(contains("bin"), rating) %>% 
  mutate(across(-rating, as.factor)) %>% 
  plot_correlation()

#this gives weird results

#now trying with rating as a factor
ingredients_df_bin %>% 
  select(contains("bin"), rating) %>% 
  mutate(across(rating, as.factor)) %>% 
  plot_correlation()

#Analysis which single ingredient is present in most recipes
df <- ingredients_df %>% 
  select(-title, -rating) %>% 
  pivot_longer(-ID, names_to = "ingredient", values_to = "value")


ing_top10 <- df %>% 
  group_by(ingredient) %>% 
  summarise(total = sum(value)) %>% 
  ungroup() %>% 
  arrange(desc(total)) %>% 
  slice(1:10)

ing_top10 %>% 
  # mutate(ingredient = fct_rev(ingredient)) %>% 
  ggplot(aes(x=reorder(ingredient, total), y=total, fill=ingredient)) +
  geom_bar(stat = "identity") +
  scale_fill_viridis(discrete = TRUE) +
  scale_x_discrete(guide = guide_axis(n.dodge=2))+
  labs(x = "Ingredient", y = "Value", title = "Total amount of recipes containing each ingredient\nTop 10")

ingredients_df_total %>% 
  select(contains("total"), rating) %>% 
  plot_correlation()

#checking some stuff about the new ingredients table
ingredients_df_total %>% 
  select(ID, title, total_ingredients) %>% 
  ggplot(aes(x=total_ingredients)) + 
  geom_bar() + geom_vline(aes(xintercept=mean(total_ingredients)),color="red", linetype="dashed", size=1)+
   scale_x_continuous(breaks = seq(1, 12, by = 1)) +
  labs(x="Number of ingredients per recipe", y = "Recipe Count", title = "Distrubution of number of ingredients per recipe")

#we notice that 117 (no NAs and RAT0, and not duplicated) recipes have 0 ingredients, let's investigate why and how that's possible
ingredients_df_total %>% 
  filter(total_ingredients==0)

#let's pick recipe ID number 1183 which should have poppy and sesame seeds according to the title
recipes %>% 
  filter(ID == 1183) %>% 
  select_if(~ any(. == 1))
#we can see that only 3 variables are equal to 1 here
recipes %>% 
  filter(ID == 365) %>% 
  select_if(~ any(. == 1))

recipes %>% 
  filter(ID == 1089) %>% 
  select_if(~ any(. == 1))

#####
#QUESTION: do we want to keep those recipes?
#####


ingredients_df_total %>% 
  filter(total_ingredients >10)

ingredients_df_full <- ingredients_df_full %>%
  filter(!total_ingredients == 0)

ingredients_df_full %>%
  filter(vegetables_bin == 1) %>% 
  select(ID, title, total_vegetables) %>% 
  ggplot(aes(x=total_vegetables)) + 
  scale_x_continuous(breaks = seq(1, 9, by = 1)) +
  geom_bar() + geom_vline(aes(xintercept=mean(total_vegetables)),color="blue", linetype="dashed", size=1)

ingredients_df_full %>%
  filter(fruits_bin == 1) %>% 
  select(ID, title, total_fruits) %>% 
  ggplot(aes(x=total_fruits)) + 
  scale_x_continuous(breaks = seq(1, 9, by = 1)) +
  geom_bar() + geom_vline(aes(xintercept=mean(total_fruits)),color="blue", linetype="dashed", size=1)

# ingredients_df_full %>% 
#   select(ID, title, total_meat) %>% 
#   ggplot(aes(x=total_meat)) + 
#   geom_bar() + geom_vline(aes(xintercept=mean(total_meat)),color="blue", linetype="dashed", size=1)+
#   labs(x="Number of meats per recipe", y = "Recipe Count", title = "Distrubution of number of meats per recipe")

#let's try to filter by recipes which contain meat to see if my functions work
ingredients_df_full %>% 
  filter(meats_bin == 1) %>% 
  select(ID, title, total_meat) %>% 
  ggplot(aes(x=total_meat)) + 
  geom_bar() + geom_vline(aes(xintercept=mean(total_meat)),color="blue", linetype="dashed", size=1)+
  labs(x="Number of meats per recipe", y = "Recipe Count", title = "Distrubution of number of meats per recipe")

#why do we still have value in 0 meats --> it was because when creating the total meat column in ingredients_df_test I did not include the general meat category